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I.  OVERVIEW 


The  title  of  the  work  reported  herein  is  quite  similar  to  that  used  for 
the  investigator's  previous  Battelle  report  [1],  Indeed,  if  the  word  "modern" 
were  omitted  they  would  be  identical.  Since  that  report  was  prepared,  Kay  & 
Marple  [2]  have  written  a  comprehensive  tutorial  paper,  "Spectrum  analysis  -  a 
modern  perspective,"  thoroughly  outlining  in  40  pages  this  topic  which  has  so 
many  potential  application  areas  to  military  radar  signal  processing,  e.g.  , 
array  processing,  clutter  suppression,  and  target  identification.  It  is 
pointed  out  in  the  literature  [3]  that  the  former  two  topics  are  basically  spa¬ 
tial  and  spectral  variations  of  the  same  problem,  namely  to  solve 

RA  =  X  (1) 

for  a  set  of  variables  (the  vector  A)  given  some  known  information  about  the 
signal  (the  autocorrelation  matrix  R)  and  some  desired  characteristics  (the 
vector  X).  Kay  &  Marple's  paper  also  includes  278  references,  many  of  which 
have  been  incorporated  into  this  report. 

The  basic  interest  in  Modern  Spectral  Analysis  (MSA)  has  centered  on  its 
superresolution  feature,  i.e.,  the  ability  to  resolve  two  closely-spaced  sinu¬ 
soids.  By  closely-spaced,  one  implies  a  frequency  separation  less  than  1/NT, 
where  N  is  the  number  of  samples  and  T  is  the  sampling  interval,  which  is  the 
accepted  resolving  capability  of  the  N-pt.  Discrete  Fourier  Transform  (DFT). 
This  claim  to  superresolution  has  been  attacked  by  Rihaczek  [4]  who  observed 
that  the  detailed  target  information  required  is  lacking  due  to  noise  or 
clutter.  In  a  rebuttal  Jackson  [5]  pointed  out  that  considerable  progress  had 
been  recently  made,  citing  among  other  things  Marple's  algorithm  [6]  which 
both -significantly  reduces  the  frequency  estimation  errors  associated  with 
Andersen's  Burg-algorithm  [7]  for  a  low  Signal-to-Noise  Ratio  (SNR)  and  vir¬ 
tually  eliminates  the  spectral-line  splitting  associated  with  a  high  SNR. 
Indeed,  in  1979  Tranter  [8]  observed  many  of  these  original  deficiencies  and 
in  a  later  study  [9]  virtually  wrote  off  this  particular  MSA  technique  because 
of  noise  problems.  This  author  in  last  years  report  made  a  similar  snap 
judgement  based  on  Tranter's  reports  which  stated  that  a  SNR  of  at  least  55  dB 
(later  reduced  by  Shumway  (10]  to  40  dB)  was  required  to  resolve  two  closely- 
spaced  sinusoids.  It  turns  out  that  the  problem  lies  in  the  technique  used  to 
estimate  the  order  of  the  prediction  filter,  and  that  it  works  for  a  SNR  as 
low  as  20  dB,  provided  the  prediction  order  is  increased.  Based  on  the 
widespread  interest  in  MSA  throughout  the  signal  processing  community,  the 
TJ.S.  Army  would  be  well  advised  to  continue  monitoring  progress  in  this  area 
for  potential  applications. 

This  report  contains  a  brief  overview  in  Section  II  emphasizing  the  Least 
Squares  Spectral  Estimation  (LSEE)  technique  associated  with  the  mc^  -Order 
Autoregressive  (AR(m))  process.  Section  III  deals  with  the  application  of  MSA 
to  the  design  of  a  Prediction  Error  Filter  (PEF)  to  suppress  clutter. 

Ironical Ly  the  presence  of  clutter  Is  one  reason  given  by  Rihaczek  for 
rejecting  MSA  techniques.  Section  IV  presents  spectral  estimation  examples 
for  AR(m)  processes  and  noisy  sinusoids  generated  from  a  computer  program 
which  is  described  in  the  Appendix.  Some  results  and  recommendations  based  on 
this  study  and  the  work  reported  in  the  literature  are  recorded  in  Section  V. 
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II.  MODERN  SPECTRAL  ESTIMATION  TECHNIQUES 


The  basic  problem  is  to  estimate  the  power  spectrum  of  a  set  of  N  samples 
(x(k)  }.  The  development  outlined  In  this  section  is  a  summary  of  the  material 
found  in  Kay  and  Marple's  tutorial  [2]  and  is  by  necessity  kept  brief.  The 
most  general  linear  model  is  the  Autoregressive  Moving  Average  (ARMA)  model 
described  by  the  linear  mc^  -order  difference  equation 


m  m 

x(n)  =  £  b(i)w(n-i)  -  £  a(k)x(n-k)  (2) 

i-0  k=l 


In  digital  filter  terminology,  a  Linear  Prediction  Filter  (LPF)  is  defined  by 
(w(n)  }  as  the  input  sequence,  (x(n)}  as  the  output  sequence,  with  {b(i)}  the 
feed-forward  coefficients  (FIR  filter),  and  (a(k)}  the  feed-back  coefficients 
(FIR  filter).  The  system  function  H(Z)  is  obtained  by  taking  the  Z-transform 
of  Eqn.  (2)  and  solving  for  the  ratio  X(Z)/H(Z). 


H(Z) 


£  b(i)Z 
i-0 


£  a(k)Z 
k-0 


(3) 


where  a(0)  A  1  and  some  b(i)  or  a(k)  can  be  zero  subject  to  the  constraint 
a(m)  #  0,  i.e.,  the  filter  is  of  order  m,  denoted  0(m).  The  power  spectral 
density  (PSD)  of  the  output  is  given  by 

1X(Z)  -  |H(Z)  |  2Pn(Z)  (4) 

2 

where  Pn(Z)  -  a  T  if  the  input  sequence  consist^  of  samples  every  Ts  from  a 
white  noise  process  of  zero  mean  and  variance  a  .  The  PSD  as  a  function  of 
frequency  is  obtained  from  Eqn.  (4)  by  replacing  Z  -*■  exp(j2xfT)  where  f  is  a 
continuous  frequency  measured  in  Hertz  (Hz). 


The  general  ARMA  model  can  be  decomposed  into  two  simpler  models,  namely 
Moving  Average  (all  a(k)  -  0,  except  a(0)  »  1)  and  Autoregressive  (all  b(i)  = 
0,  except  b(0)  =  1)  which  by  Eqn.  (3)  can  also  be  termed  the  all-zero  and  all¬ 
pole  models,  respectively.  The  autoregressive  (AR)  model  is  the  subject  of 
primary  interest  in  this  study.  It  follows  from  Eqn.  (4),  assuming  white 
noise  input,  that  the  AR  PSD  is 


P(f)  - 


02T 


1  +  £  a(k)  exp(-j2ufkT) 

k- 1 


(5) 


The  numerator  is  a  scale  factor;  consequently,  the  relative  PSD  is  completely 
determined  by  knowledge  of  the  parameter  set  (am(k)>  of  the  AR  process  of 
order  m,  AR(m) ,  with  the  understandings  that  the  parameter  subscript  m  can  be 
dropped  whenever  obvious  and  am(0)  -  1  for  all  m.  Furthermore,  the  final 
value  in  each  set  is  known  in  the  literature  as  the  reflection  coefficient 
( K^) >  i .e«  , 


2 


=  Vm> 


(6) 


The  FIR  digital  filter  with  weights  (a^k)},  x(n)  as  input  and  w(n)  as  output 
is  called  a  Prediction  Error  Filter  (PEF).  The  obvious  question  remains, 
given  a  sample  set  (x(k)}  how  does  one  obtain  {a(k)}  in  order  to  determine 
P( f )  from  Eqn.  ( 5) . 

A.  Yule-Walker  Equations  (YWE) 

The  original  developments  for  finding  the  (a(k)}  for  the  Mth  -order 
LPF  assumed  knowledge  of  the  autocorrelation  values,  which  in  theory  are  given 
by 


R(k)  *  E{x(u+k)x*(n) }  (7) 

where  E(  }  implies  statistical  expectation  and  the  superscript  (*)  complex 
conjugation.  The  YWE 


R(k) 


M 

E  a(i)R(k-i)  k  >  0 

i=l 


2 

E  a(i)R(-i)  +  a  k  -  0 
i-1 


(8) 


follow  by  substituting  Eqn.  (2^  into  Eqn.  (7)  and  recognizing  that 
E  (n(n+k)x*(n)  }  is  non  zero  (  =  o  )  only  for  k=0.  The  AR  parameter  set  is 
obtained  by  selecting  M  equations  from  the  semi-infinite  set  of  Eqn.  (2)  and 
solving  for  a  from  the  k=0  equation.  In  matrix  form  this  is 


[R(i-k) ] (M,M)  x  [a(i)](M>1)  =  -  [R(i)](M,i)  (9) 

where  subscript  (x,y)  implies  a  matrix  of  x  rows  and  y  columns.  The  matrix 
[R(i-k)]  is  Herml tian  and  Toeplitz  i.e.,  R(i,k)  =  R(i-k)  «  R*  (k-i). 
Rearranging  Eqn.  (9)  and  including  the  row  for  a  yields 


[R(i-k) l(M+i ,M+i)  (I  a(l)  ...  a(M) ] '(M+1 


[ a  0  0  ...  0] 


(M+1,1) 


(10) 


where  prime  (  ')  implies  a  vector  transpose.  The  YWE  given  by  Eqn.  (10)  can  be 
solved  for  fa(i))  by  inverting  the  matrix  [R(i-j)]  provided  the  first  (M+l)- 
lags  of  the  autocorrelation  are  known.  Because  of  the  Toeplitz  structure,  the 
standard  matrix  inversion,  requiring  0(M  )  operations ,^can  be  replaced  by  the 
Levinson-Durbin  Equation  (LDE),  which  only  require  0(M  )  operations.  The 
technique  is  to  recursively  solve  the  AR(m)  parameter  set  (am(k)}  for  m  = 
2,3,...,M  using  the  algorithm 
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(11) 


_  p  m- 1 

am(m)  A  Kjj,  =  cr^Wm)  +  T.  a^.^i)  R(m-i)] 


am(i)  =  am_ ! ( i )  +  a^.jCm-i) 


^  =  (1  -  K  |2]  a2^ 


subject  to  the  Initial  conditions 


a2  =  [1  -  |  Kx  |  2]  R(0) . 

Since  the  desired  order  (M)  of  the  LPF  is  not  known  a  priori,  Eqn.  (11)  can  he 
solved  for  successively  higher-order  values  of  tn  until  the  modeling  error 
om  is  reduced  to  an  a^ceptab^e  value.  Further,  since  |  Km  |<  1,  it  follows 
from  Eqn.  (11)  that  om4.^  <  om  ;  moreover,  if  the  (x(k)}  are  truly  samples  of 
an  AR(M)  process,  the  modeling  error  reaches  an  irreducible  minimum  o m  .  The 
poles  of  A(Z),  the  denominator  of  Eqn.  (3),  can  be  shown  to  lie  within  the 
unit  circle,  i.e.,  the  LPF  is  stable,  or  marginally  stable  (poles  on  the  unit 
circle)  if  the  AR  process  consists  only  of  sinusoids. 

B.  Maximum  Entropy  Spectral  Estimation  (MESE) 

The  concept  of  AR(m)  modeling  and  associated  LPF  design  has  been  shown 
to  be  equivalent  to  the  Maximum  Entropy  Method  (MEM)  propost d  originally  by 
Burg  [11,  12],  provided  the  autocorrelation  lags  are  known  and  the  random  pro¬ 
cess  is  Gaussian.  Burg  suggested  that  extrapolation  of  the  m  known  correla¬ 
tion  lags  should  be  performed  in  such  a  manner  that  the  corresponding  time 
series  has  maximum  entropy,  i.e.,  the  PSD  should  be  as  white  as  possible  given 
the  known  lags.  The  MESE  is  given  by  Eqn.  (5)  with  the  (a(k)}  obtained  from 
solution  of  Eqn.  (10). 

C.  Estimating  the  AR  Parameters 

Unfortunately,  in  a  radar  system  the  { R( k ) }  values  are  unknown  and 
all  one  can  begin  with  are  the  N  time  samples  (x(k)>  from  an  assumed  AR(m) 
process.  Although  one  might  estimate  the  desired  correlation  samples  by  the 
biased  estimator. 


R(k)  =  _2_  T.  x(n)x  (n+k)  (13! 

N  n=1 

it  has  been  shown  for  short  data  records  (N  100)  the  YWE  produces  poor 
spectral  estimates. 


FREQ.  0.11250000E+00  MAX.  0.0000000 E *0 


Figure  4.  AR(4)  spectral  estimate  using  MLSA. 

B.  Noisy  Sinusoid  Estimation 

The  first  problem  to  be  examined  is  a  noisy  sinusoid  given  by 

y(k)  =  A  cos(  ek  +  4>)  +  n(k)  (41) 

where  A  2  nf CT  and^nCk)  is  a  sample  from  a  complex  Gaussian  process  with  zero 
mean  and  variance  a  in  either  the  real  or  imaginary  component.  The  SNR  (not 
In  dB)  is  given  by 

2 ,  2 

a  =  A  /2o  (A2) 

and  Is  applicable  for  both  real  and  complex  data.  Swingler  showed  [45]  that 
the  signal  portion  of  Eqn.  (41)  when  estimated  by  the  ABA  has  a  frequency 
estimate  which  differs  from  fr  bv 


T 


Figure  3.  AR(2)  spectral  estimates  using  ABA. 

The  second  example  is  a  40  sample  AR(4)  process  with  coefficients 
f2.7607,  -3.8106,  2.6534,  -0.9238}  and  CNR  =  50  dB.  the  spectral  estimate 
which  results  from  setting  the  maximum  allowed  order  (MM)  equal  to  four  in  the 
MLSA  is  shown  in  Figure  4.  This  spectrum  is  suggestive  of  that  produced  by 
broadband-type  (bandlimited  white  noise)  jamming  and  presumably  could  be 
suppressed  in  a  fashion  similar  to  clutter  using  a  PEF.  The  algorithm  ter¬ 
minated  at  M=4  because  the  ratio  of  f inal-to-initial  prediction  error  energies 
was  below  T0L1  =  0.001.  Essentially,  an  identical  result  for  M=4  was  obtained 
with  the  ABA;  however,  the  three  PEF-order  indicators  tabulated  in  Figure  A-2 
all  predicted  M  >  4.  FJlrych  &  Bishop  observed  that  the  histogram  of  absolute 
minimum  FPE-order  indicators  obtained  from  100  estimates  of  AR(2)  and  AR(4) 
processes  was  rather  broad,  but  improved  with  significiant  numbers  at  n=M  if 
MM  was  restricted  to  N/3.  They  also  determined  that  the  histogram  of  first 
FPE  minimum  (vs.  absolute  minimum)  was  quite  similar,  an  approach  to  selecting 
M  which  is  commonly  used. 
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A.  AR  Model  Estimations 


Two  AR  models  specified  by  Ulrych  &  Bishop  [15]  for  real  data  were 
used  to  verify  the  two  algorithms  and  demonstrate  the  applicability  of  this 
technique  to  wideband  spectral  estimation.  The  first  example  was  a  20  sample 
AR( 2)  process  with  coefficients  {0.75,  -0.50}  and  CNR  =  50  dB  producing  a 
theoretical  spectrum  which  looks  somewhat  like  weather  clutter.  The  results 
obained  with  MLSA  for  M=3  and  M=4  are  shown  in  Figure  2.  The  M=3  solution  is 
obviously  not  a  good  estimate;  however,  the  M=4  solution  is  essentially  iden¬ 
tical  to  that  shown  in  Ulrych  &  Bishop.  Results  obtained  from  the  ABA  for 
M=3,  4  and  11  are  shown  in  Figure  3.  Obviously,  the  ABA  does  not  produce  as 
accurate  a  prediction.  Moreover,  the  M=ll  example  shows  the  effect  of 
overestimating  M  excessively.  A  change  of  N  from  20  to  50  samples  had  no 
noticable  effect  on  the  spectral  estimates. 


Figure  2.  AR(2)  spectral  estimates  using  MLSA. 
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Gibson  and  Haykin  observed  that  the  number  of  sections  required  is 
at  least  equal  to  the  number  of  spectral  peaks,  equal  for  sharp  peaks,  but 
typically  two  or  three  for  each  broad  peak  such  as  clutter  produces.  They 
demonstrated  the  adequacy  of  a  notch  filter  (M=l)  for  a  sinusoid  and  the  effect 
of  changing  M  for  a  Gaussian  shaped  clutter  spectrum,  see  Eqn.  (27). 

Increasing  M  from  1  to  3  essentially  broadened  the  stopband  to  match  the 
clutter  bandwidth  while  additional  stages  did  little  except  to  flatten  the 
passband  gain.  They  also  illustrated  how  to  select  tj  when  the  clutter  extent 
is  tc  »  100  and  signal  duration  ts  =  20.  Initial  clutter  modeling  showed  that 
decoupling  occurred  when  tj  =  25  for  W  =  0.9.  It  follows  from  Eqn.  (39)  that 
a  =  2.64,  and  from  their  second  choice  of  W  =  0.95  with  signal  and  cluster 
present,  that  the  new  tj  =  51  is  consistent  with  setting  tj  =  tc/M  where  a 
second-order  (M*2)  PEF  was  specified.  A  flow  chart  for  implementating  the 
recursive  reflection  coefficient  computations  for  either  method  is  provided 
[43]. 


In  two  conference  papers  Gibson,  e_t.  al_.  ,  studied  the  adaptive  lat¬ 
tice  PEF  using  simulated  [39]  and  actual  [44]  radar  clutter  data  in  complex 
form.  The  earlier  paper  (1979)  with  simulated  data  employed  Method  I  and 
defined  W  =  0.95  and  M  *  2  based  upon  a  large  number  of  test  cases.  It  was 
observed  that  W  is  a  function  of  the  number  of  clutter  samples  (N)  with  W  = 
0.98  for  N  =  250  vs.  0.95  for  N  =  100.  Published  filter  responses  were  shown 
to  be  a  function  of  CNR  with  little  variation  for  CNR  _>  10  dB.  When  expressed 
in  terms  of  SIR,  the  average  enhancement  exceeded  20  dB  for  ’’average  storm" 
conditions  and  16  dB  for  "serve  storm”  conditions.  The  latter  reflects  a 
situation  where  the  conventional  MTI  would  not  properly  suppress  the  clutter 
due  to  its  increased  bandwidth.  The  later  paper  (1981)  tested  the  performance 
of  Che  lattice  PEF  employing  Method  II  using  actual  radar  returns  from  sta¬ 
tionary  and  nonstationary  clutter  with  and  without  targets  present.  The 
improvement  ratio  (SCR-out  over  SCR-in  for  MTI  or  PEF)  was  always  positive  for 
the  PEF,  ranging  from  2  to  8  dB,  whereas  the  range  was  -17  to  +16  dB  for  the 
MTI  with  poor  performance  associated  with  adverse  weather  conditions.  The 
radar  used  employed  pulse  staggering  with  four  unique  interpulse  spacings.  In 
order  to  employ  this  data  with  the  lattice  filter,  the  signal  returns  were 
subdivided  into  four  groups  each  having  five  target  samples  and  25  clutter 
returns.  The  subsets  were  filtered  independently  and  the  displayed  graphs 
represented  the  average  of  the  four  subsets  for  input,  MTI  and  lattice  PEF 
outputs.  The  choice  of  adaptive  weighting  factor  (U=0.6)  and  filter  order 
(M=5)  were  determined  only  after  extensive  testing,  probably  required  by  the 
small  number  of  samples  (N=20)  available. 

IV.  SPECTRAL  ESTIMATOR  PERFORMANCE 

The  maximum  entropy  and  least  squares  spectral  estimators  were  described 
in  Section  II  and  their  applicability  to  clutter  suppression  in  Section  III. 
The  purpose  of  this  chapter  is  to  present  details  regarding  their  applicabi¬ 
lity  to  the  spectral  estimation  of  AR(m)  processes  and  noisy  sinusoids,  both 
single  component  and  two  sinusoids  spaced  closer  than  the  DFT  resolution 
width.  The  computer  algorithms  ABA  and  MLSA  are  used  for  the  MESE  and  LSSE , 
respectively.  A  computer  program  which  can  generate  noisy  complex  data 
samples  and  estimate  the  spectrum  using  either  algorithm  is  described  in  the 
Appendix  along  with  some  typical  examples. 
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exponentially  weighted  and  eventually  eliminated.  Method  II  updates  the 
crosspower  (numerator)  and  autopower  (denominator)  expectations  separately 
which  yields 


where 


*W»> 


vm+l(n) 

ym+i(n> 


(36) 


vm+l<n) 

ym+^n)  ,  U 


U  •  vm+1(n-l)  -  2fm(n)  .  bm(n-l) 
ym+l<n_l)  +  |fm<n>  |  2  +  IV11-1)  |  2 


(37) 


with  U  being  the  new  adaptive  weighting  factor  and  vm+j(0)  **  ym+i  (0)  “  0. 

The  forward-backward  prediction  errors  were  defined  in  Eqn.  (15).  Method  II 
has  the  advantage  of  not  assuming  constant  power  and  should  work  best  with 
nonstationary  signals.  However,  it  is  more  complicated  computationally  and 
has  a  more  complex  convergence  relationship.  They  also  showed  that  Methods  I 
and  II  were  special  cases  of  Griffiths’  algorithm.  It  was  further 
demonstrated  that  for  Method  II  Kj  converged  to  its  optimal  value  more  rapidly 
with  a  corresponding  faster  decrease  in  |f,(n)  |.  The  problem  of  coefficient 
decoupling  was  also  addressed  and  shown  to  &e  a  significant  consideration. 

Decoupling  is  a  problem  because  initially  all  {K^}  in  the  Mth  -order 
PEF  adapt  to  the  strongest  spectral  component.  Eventually  Kj  adapts  properly 
and  the  remaining  M-l  coefficients  attack  the  next  strongest  component,  etc. 
Furthermore,  residual  coupling  remains  after  the  major  decoupling  occurs,  thus 
interfering  to  some  extent  with  convergence  properties  of  later  stages.  The 
first  stage  decoupling  time  (tj)  is  proportional  to  n,  i.e. , 


tj  »  an  (38) 

where  a  is  a  dimensionless  constant  associated  with  a  specific  signal  type. 

For  method  I,  it  follows  from  Eqn.  (35)  that 

W  »  exp(-a/tj)  .  (39) 

Unfortunately,  a  similar  expression  does  not  exist  for  the  weighting  factor  U 
of  Method  II.  Assuming  the  incoming  signal  consists  of  both  clutter  with 
duration  tc  and  signal  with  considerably  smaller  duration  tg,  then 

tg  <  ci  <  tc  •  (40) 

As  the  filter  order  increases  tj  is  set  closer  to  tg,  and  for  very  large  tc, 
the  upper  limit  is  typically  divided  by  an  integer  ranging  from  one  to  M,  thus 
insuring  time  for  higher-order  stages  to  converge. 
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rejected  in  favor  of  the  original  ABA  because  it  required  three  times  as  many 
complex  multiplications  for  this  application.  He  also  rejected  Andersen’s 
suggested  modification  [13],  i.e.  ,  Eqn.  (18)  on  the  basis  of  sensitivity  to 
computational  errors.  The  current  RADC  system  can  process  128  range  cells  and 
allows  32  updates  of  the  weight-sector  components  during  this  dwell  time. 
Although  currently  operating  with  the  Burg  algorithm,  the  adapative  processor 
is  reprogrammable.  Documentation  of  the  RADC  system  should  be  available 
shortly. 

C.  Adaptive  Lattice  PEF 

Gibson  and  Haykin  in  a  series  of  paper  [39]  show  that  clutter 
suppression  can  be  achieved  using  an  adaptive  PEF  with  the  lattice  structure 
as  opposed  to  the  more  conventional  Tapped  Delay  Line  (TDL)  used  to  implement 
Eqn.  (32).  The  lattice  PEF  structure  was  shown  in  Figure  1  and  the  recursion 
relations  were  given  by  Eqn.  (15).  Burg's  (harmonic-mean)  algorithm  states 
that  the  optimum  values  of  the  reflection  coefficients  {K^,}  are  given  by  the 
ratio  of  the  expectations  of  the  negative  crosspower  to  the  mean  value  of  the 
two  prediction  error  autopowers  for  that  PEF  stage  (m).  Makhoul  [40]  descri¬ 
bes  a  variety  of  alternatives  to  the  harmonic  mean  definition  for  calculating 
{K,,,},  but  concludes  that  Burg's  is  preferred  "because  it  minimizes  a  reaso¬ 
nable  and  well-defined  error  criterion",  a  view  supported  by  Gibson  [41].  If 
adaptation  were  not  required,  i.e.,  the  input  signal  was  from  a  stationary 
process,  the  {K^,}  for  Burg's  algorithm  are  given  by  Eqn.  (17). 

When  clutter  is  to  be  suppressed  in  either  the  presence  or  absence 
of  a  doppler-return  signal,  it  follows  that  the  PEF  coefficients  must  be 
updated  routinely,  i.e.,  adapted  to  changes  in  the  radar  return.  Griffiths 
[42]  compares  the  lattice  and  TDL  filters  from  the  viewpoint  of  adaptive 
filtering  and  concludes  that  the  lattice  enjoys  convergence  advantages  at  the 
expense  of  increased  computations.  The  lattice  structure  for  stationary  pro¬ 
cesses  permits  the  independent  optimization  of  each  which  carries  over  to 
some  degree  for  the  adaptive  case.  It  is  also  possible  to  trade  off  adap¬ 
tation  rate  and  prediction  error  with  filter  length.  Gibson  and  Haykin  [43] 
studied  the  properties  of  adaptive  lattice  PEF  and  proposed  two  methods  for 
recursive  estimation  of  1^  on  an  incoming  sample-by-sample  basis.  Method  I 
adds  a  correction  term  to  the  previous  value  which  yields 

K^^n)  =  W  .  l^+^n-l)  +  Atn(n)  .  fm(n)  .  bm(n-l),  (34) 


where  fm(n)  and  bm(n)  are  given  by  Eqn.  (15),  and 


with 


Am(n)  - 


-2( 1-W) _ 

+  |Vn-l) 


W  =  exp(-l/n)  , 


(35) 


(36) 


where  n  is  the  length  (in  samples)  required  for  the  filter  to  adapt.  The 
adaptive  weighting  factor  (W)  is  designed  to  insure  that  old  estimates  are 


Incorporating  the  adaptive  whitening  fiiter  and  the  DFB  in  one  step 
has  been  reported  by  Sawyers  [38,  25],  The  resulting  coefficients  are 
obtained  from  Eqn.  (31)  with  R“^  obtained  from  Eqn.  (25)  using  the  MEM 
algorithm.  The  coefficient  set  obtained  from  N  samples  in  one  dwell  is  used 
to  define  the  adaptive  DFB  coefficients  for  the  next  dwell  in  the  same  range 
celi.  Should  a  target  be  detected  in  that  cell,  the  update  algorithm  is  inhi¬ 
bited  so  as  not  to  suppress  the  target  spectrum.  Sawyers  presents  performance 
results  for  a  pulse-doppler  radar  (N  =  32,  M  =  15)  using  a  simulated  inter¬ 
ference  spectrum  including  ground,  rain  and  broadband  jamming  with  clutter-to- 
noise  ratios  of  50,  20,  and  30  dB ,  respectively.  The  simulation  assumes  unit 
signal  voltage  and  unit  thermal  noise  power  per  returned  pulse.  Optimally,  a 
SIR  of  14.8  dB  is  possible  for  a  particular  doppler  fiiter  (fnT  »  0.375)  with 
the  actual  design  indicating  14.0  dB.  The  improvement  is  definitely  a  func¬ 
tion  of  Doppler,  as  a  second  Doppler  choice  (fnT  -  0.5625)  just  below  the 
broad-band  jamming  has  an  optimal  SIR  of  11.9  dB  with  the  actual  result  some 
2.5  dB  less.  It  is  suggested  that  additional  dwell  averaging  be  invoked  in 
such  cases,  i.e.,  the  reflection  coefficient  Km  given  by  Eqns.  (17)  and  (18) 
becomes  for  D-dwell  averaging 


D 

£  Num( i ) 


£  Den(i) 
i-1 

where  Den(i)  is  defined  by  Eqn.  (18)  and  Num(i)  is  the  summation  expression  in 
Eqn.  (17).  Only  limited  results  were  presented  for  actual  radar  data  and  no 
attempt  was  made  to  determine  the  SIR. 

Nitzberg  [23]  reported  an  adaptive  DFB,  similar  to  Sawyers,  designed 
by  General  Electric  for  the  USAF  Rome  Air  Development  Center,  (RADC) .  This 
system  was  implemented  with  the  Burg  algorithm  and  did  a  MESE  for  N  *  16.  It 
appears  that  the  estimated  clutter  spectrum  for  one  range  cell  was  used  to 
filter  the  clutter  in  the  next  cell  as  opposed  to  the  same  cell  one  dwell-tirae 
later.  There  was  no  discussion  of  dwell  averaging  to  improve  the  estimated 
spectrum.  Nitzberg  did  study  the  effects  of  a  ’’Mismatch  Loss"  due  to  the 
AR(M)  spectrum  not  matching  the  true  spectrum  by  investigating  the 
Clutter-to-Noise  Ratio  (CNR)  after  filtering  vs.  filter-order  (M)  for  various 
input  CNR.  The  simulated  clutter  had  a  Gaussian  spectrum,  from  which  the 
correlation  matrix  can  be  determined,  hence  the  optimum  weights  from  Eqn.  (31) 
and  residual  CNR  decreased  to  within  a  few  dB  of  the  theoretical  optimum.  Mow 
rapidly  it  converged  to  the  optimum  depended  on  the  original  CNR,  e.g.  ,  if  the 
input  CNR  is  35  dB  the  residual  is  4  dB ,  whereas  for  65  dB  the  residual  is 
roughly  14  dB,  which  is  some  10  dB  above  the  optimum  weight  result.  The 
mismatch  loss  is  also  shown  to  be  dependent  upon  the  normalized  clutter  stan¬ 
dard  deviation  ( a  T)  and  the  relative  separation  of  the  clutter  mean  from  the 
doppler  filter  center  frequency,  increasing  with  the  former  and  decreasing 
with  the  latter  as  they  increase  in  value.  Nitzberg  also  discussed  an 
"Estimation  Loss”  due  to  the  AR(M)  parameters  being  estimated  from  data 
samples  as  opposed  to  known  values  of  the  R-matrix.  He  showed  that  the  esti¬ 
mation  loss  tends  to  increase  with  the  LPF  order  thus  making  it  difficult  to 
determine  the  optimum  order.  The  more  recent  Marple  algorithm  [6]  was 


B.  Adaptive  Doppler  Signal  Processor 

Other  researchers  [3,  34,  35]  had  proposed  designing  filters  to 
maximize  the  Signal-to-Clutter  Ratio  (SCR)  as  opposed  to  the  conventional 
match  filter,  which  maximizes  the  SNR,  where  the  noise  is  assumed  to  be  white 
Gaussian.  The  problem  was  that  the  clutter  characteristics  had  to  be  spe¬ 
cified  ji  prloi .  Two  examples  where  this  philosophy  was  employed  were  the 
Hughes’  Aircraft  Co.  designs  of  a  Doppler  Filter  Bank  (DFB)  for  the  TPQ-37 
artillery  locator  [36]  and  the  U.  S.  Army  Missile  Command's  Quiet  Radar. 
Basically,  a  conventional  matched  filter  can  be  implemented  for  a  given 
Doppler  shift  (fn)  using  a  Discrete  Fourier  Transform  (DFT),  whereas  a  fiiter 
matched  to  a  specified  Signal-to-Interference  Ratio  (SIR)  is  obtained  by  maxi¬ 
mizing 


SIR  ,  V  ¥  •  <2,> 

wt  Rw 

Assuming  uniform  transmit  weights,  s  is  s  steering  vector  defined  by 

£(n)+  =  [exp(-j  ir(N-l)fnT)  ,  exp(-j  x(N-3)fnT . exp  (+j*(N-l)fnT)]  ,  (30) 

and  R  is  the  N  x  N  matrix  of  correlation  samples.  The  vector  of  filter  coef¬ 
ficients  w  can  be  obtained  by  solving 

w(n)  -  R-1£(n)  ,  (31) 

where  the  ”n"  is  identified  with  the  Doppler  frequency  fn.  The  solution  of 
Eqn.  (31)  depends  on  finding  R-1 ,  which  was  discussed  previously  in  Section 
II. 1,.  Such  an  Inverse  technique  has  a  computational  complexity  of  (|(N)  to 
0(N  )  as  opposed  to  direct  inversion  by  Gaussian  elimination  of  0(N  ). 

Another  technique  for  estimating  R”1  involves  Kalman  filtering  which  has  a 
computational  complexity  0(N  ).  This  approach  was  used  by  Bowyer  ££  al_  [37] 
in  designing  an  adaptive  clutter  filter  for  a  problem  involving  ballistic 
missiles  where  booster-tank  break  up  upon  atmospheric  reentry  caused  a  false 
target.  They  observed  that  their  Kalman  estimator  approach  was  far  superior 
to  the  earlier  FIR  filters  which  used  feedback  to  adjust  the  weights  and  hand 
a  long  convergence  time  on  the  order  of  NT  (20  <  N  <  30).  Their  technique 
assumed  an  AR(M)  clutter  model  (2  £  M  £  4)  with  the  spectral  coefficients 
(a(i)}  obtained  from  the  Kalman  filter  used  to  construct  a  whitening  filter 


HW(Z)  «  1  -  T.  a(i)  Z-1  , 

i-1  (32) 

whose  structure  can  be  readily  changed  with  the  clutter.  The  entire  N  samples 
were  then  whitened  and  sent  to  a  conventional  matched  fiiter  with  coefficients 
adjusted  to  reflect  signal  whitening.  No  attempt  was  made  to  integrate  the 
whitening  process  with  the  DFB.  It  was  also  pointed  out  that  the  estimator 
memory  is  a  drawback  to  using  the  Kalman  filter  algorithm.  Consequently,  if 
the  clutter  is  not  extended  in  time,  the  estimates  do  not  change  rapidly  and 
It  becomes  necessary  to  introduce  some  mechanism  such  as  exponentially  fading 
memory  to  insure  proper  adjustment. 
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which  can  be  approximated  [27]  with  an  error  |  e  )  <  2.7  x  10“  3  by 

Pg( F)  =  [C(0)  +  C(2)F2  +  C(4)F4  +  C(6)F6]“1  ,  (28) 

with  (C(n)  }  -  (2.490895,  1.466003,  -0.024393,  0.178257).  It  should  be  noted 
that  F  represents  the  actual  frequency  normalized  with  respect  to  a  parameter 
which  depends  both  upon  the  transmit  frequency  and  the  nature  of  the  clutter. 
The  important  point  is  that  both  forms  can  be  considered  as  all-pole  functions 
and  hence  modeled  as  AR(M)  processes  where  M  is  n/2  for  the  Butterworth  model 
(typically  2  <_  n  4)  and  three  (3)  for  the  Gaussian-approximation  model  of 
Eqn.  (28).  Since  n  in  Eqn.  (26)  is  not  required  to  be  an  integer  and  the 
Gaussian  model  is  an  approximation,  a  realistic  M-value  would  be  somewhat 
larger.  Obviously,  larger  M-values  imply  more  computation  and  peaky  spectral 
approximations. 

Hawkes  and  Haykin  [28]  developed  a  computer  model  for  I&0  channel 
(complex-data)  clutter  samples  which  assumes  a  large  number  of  elemental 
dipole  scatters  arranged  in  a  two-dimensional  (range  and  azimuth)  array.  The 
model  Included  Probability  Density  Functions  (PDFs)  to  describe  the  dipole 
rotation  frequency,  doppler  shift  due  to  clutter,  and  a  factor  to  describe  the 
antenna  pattern. 

A.  Clutter  Spectral  Estimation 

Kesler  and  Haykin  [29,  30]  used  Burg's  MEM  technique  to  estimate  the 
clutter  PSD.  The  aforementioned  computer  data  was  used  assuming  a  steady 
dipole  frequency  and  describing  the  doppler  variations  by  a  Gaussian  PDF  of 
zero-mean  and  standard  deviation  a.  Akaike’s  FPE  algorithm  was  used  to  deter¬ 
mine  the  appropriate  choice  for  M.  They  also  generated  samples  of  AR(2)  and 
AR( 4)  processes  (real-valued)  which  yielded  rather  poor  MESE,  looking  neither 
like  the  theoretical  curves  or  each  other  as  M  was  varied  above  and  below  the 
FPE  estimate.  However,  when  complex  data  was  used,  either  from  the  computer 
model  or  an  AR  process,  the  resulting  spectra  were  relatively  Insensitive  to 
the  choice  of  M  or  N.  It  should  be  observed  that  N  was  typically  quite  large 
(  64)  for  the  AR  process  and  that  the  results  were  sensitive  to  start-up 

effects,  i.e.,  the  number  of  AR-model  samples  ignored.  However,  the  computer- 
model  data  results  looked  similar  to  the  theoretical  spectrum  for  N  as  low  as 
16.  By  comparison  spectra  obtained  by  applying  Welch's  periodogram  method 
(Fourier  analysis)  to  the  same  data  yielded  unacceptably  broad  spectra  for  N 
64.  In  later  papers  [31,  32]  the  MEM  approach  was  applied  to  true  radar 
clutter  (complex  data)  including  ground,  weather  and  birds.  The  ground  and 
weather  clutter  spectra  were  quite  insensitive  to  N  provided  it  was  In  the 
range  10  to  30  samples  or  greater.  The  bird-spectra  were  quite  variable  with 
N  in  agreement  with  variable  spatial  and  temporal  distributions  within  flocks. 
Their  conclusions  were  that  good  MESE  could  be  obtained  for  a  paricular  range 
bin  using  only  one  look.  They  also  suggested  that  MEM  might  be  used  for  on¬ 
line  classification  of  various  clutter  sources.  Hawkes  and  Haykin  in  an 
earlier  paper  [33]  had  proposed  an  adaptive  Moving  Target  Indicator  (MTI)  the 
coefficients  of  which  were  chosen  based  on  a  decision  as  to  what  type  of 
clutter  was  encountered.  This  decision  was  based  on  applying  the  autocorrela¬ 
tion  of  the  incoming  signal  to  a  decision  algorithm  and  comparing  the  output 
against  a  set  of  precomputed  clutter  characteristics. 
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where  (f)  implies  a  complex-conjugate  transpose  operation.  The  inverse  matrix 
is  both  Herraitlan  (R~l  =  (R”l)t  )  and  persymmetric  (symmetric  about  both 
diagonal  directions).  The  matrix  A  is  N  x  N  with  the  first  N  -  M  columns  con¬ 
sisting  of  (aM(k)}  with  {0,1,2,  ...,(N-M)  leading  zeros  and  (N-M-l),  ...,2,1} 
trailing  zeros  respectively,  plus  M  columns  consisting  of  (a^Ck)},  m  =  0  *  M 
arranged  with  the  proper  number  of  leading  zeros  (N-M  for  ra=M  up  to  N-l  for 
m=0).  The  matrix  D  is  an  N  x  N  diagonal  matrix  with  1 / EM  for  the  first  N-M 
diagonal  elements  and  {l/E^_^ , . . . ,  1/Ej,  1/Eg}  for  the  remaining  M  elements. 
Details  regarding  the  derivation  of  Eqn.  (24)  beginning  with  the  YWE  are  also 
found  in  Herring  [22],  Nitzberg  [23]  observed  that  it  should  be  possible  to 
compute  R”l  using  only  the  PEF  set  for  m=M  and  suggested  a  technique  for 
complex-data  based  on  Siddiqui  [24],  Although  Nitzberg's  algorithm  was  not 
published.  Sawyers  [25]  independently  presented  such  a  solution  assuming  M  <_ 
N/2  and  showed  that  {Z( i , j ) }  the  elements  of  R"1  are  given  by 

Z(i,j)  =  a( i-l)a*( j-1)  +  Z( i— 1 ,j-l)  (25a) 

where  i  =  1,2,  ...,  min  (N/2,M+1)  and  j  =  i,i+l,  ...,  M+l;  or 

Z(i  ,  j )  -  Z(i-l,j-l)  (25b) 

for  i  =  2,3,  ...,  rain(N/2,M+l)  and  j  *  M+2,M+3,  ...,  provided  N  2(M+2) ,  also 
for  i  =  M+2,M+3,  ...,  N/2  and  j  =  i,  i+1 ,  ...»  i+M;  or 

Z(i,j)  =  0  (25c) 

for  i  =  1,2,  ...,  N/2  and  j  *  i+l+M,  i+2+M,  ...»  N;  or 

Z(j,i)  =  Z*(i  ,  j )  ( 25d) 

for  i  =*  1,2,  ...,  N/2-1  and  j  =*  i+l,i+2,  ...»  N/2.  The  elements  of  the 
lowerhalf  of  R“1  could  be  obtained  from  the  fact  that  R~*  is  Hermitian  and 
persymmetric.  In  reality,  the  upper-half  of  R~1  is  sufficient  for  the  clutter 
application  described  in  Section  III. 

III.  CLUTTER  SUPPRESSION 

Clutter  is  a  general  terra  describing  radar  returns  from  objects  other 
than  targets  of  interest  e.g.  ,  ground,  rain  clouds,  or  bird  reflections.  It 
does  not  Include  energy  received  from  deliberate  attempts  to  interfere  with 
radar  performance  such  as  chaff  or  active  jamming.  Zehner  and  Currie  in  a 
recent  study  [26]  suggest  that  the  best  PSD  models  are  Butterworth 

Pb(F)  =.  _  (26) 

1+Fn 

and  Gaussian 

Pk(F)  =  J_  exp(-F2/2)  , 

>/z7  (27) 
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Inspection  of  Eqns.  (22)  and  (23)  show  that  a  penalty  is  paid  for  selecting 
m  >  M,  since  decreases  slower  than  the  compensating  factor  or  term 
increases.  Ulrych  and  Bishop  [16]  have  studied  the  Akaike  FPE  criterion  for 
both  YWE  and  ABA.  Histograms  of  the  minimum  FPE  location  for  100  realizations 
of  AR(2),  (N=20)  and  AR(4),  (N*4Q)  data  sets  indicated  peaks  when  m»M  and  m=N. 
When  the  order  was  constrained  (m  N/3)  a  single  peak  near  m*M  occurred. 

Their  studies,  supported  later  by  Jones  [17],  indicated  that  with  N/2  the 
first  minimum  in  Eqn.  (22)  should  be  used  to  determine  the  LPF  order  for  an 
AR(M)  process.  If  the  data  comes  from  noisy  sinusoids  a  clear  minimum  for  the 
FPE  does  not  exist,  and  the  first  local  minimum  could  easily  produce  a  poor 
spectral  estimate.  Landers  and  Lacoss  [18]  studies  all  three  order-predictors 
which  showed  that  the  CAT  technique  produced  a  sharp  minimum  for  noise-free 
sinusoids  when  compared  to  the  Akaike  methods.  However,  there  was  little 
distinction  among  the  three  for  the  same  sinusoid  plus  50%  uniform  noise.  It 
did  appear  that  a  single  minima  occurred  in  predictor-algorithm  vs.  order- 
selected  plots  when  processing  a  complex  sinusoid  as  contrasted  with  either 
the  real  or  imaginary  components.  Unfortunately,  it  is  often  the  size  and 
location  of  narrowband  spectral  components  which  are  of  interest,  whereas,  the 
order-predictor  algorithms  are  based  upon  the  entire  spectrum  including  any 
noise  present. 


H.  Compensating  Observation  Nose 

When  observation  noise  is  added  to  samples  of  an  AR  process,  the 
result  is  realiy  an  ARMA  process  and  spectral  estimates  derived  from  LSSE  and 
MESE  are  inaccurate,  becoming  worse  as  the  Signal-to-Noise  Ratio  (SNR) 
decreases.  Unfortunately,  ARMA  modeling  techniques  are  not  well  developed, 
particularly  for  short  data  records  [2,  p.  1397],  A  second  approach  is  to  use 
the  AR  modeling  but  let- m  =  N/3.  Kay  [19]  has  shown  that  it  is  possible  to 
detect  a  noisy  sinusoid  (SNR  =  OdB)  with  m=32  (N-100).  In  a  subseqent  paper 
[20]  Kay  investigated  a  second  approach,  namely  compensating  the  R-matrix 
directly  by  subtracting  00^1,  where  a  is  a  decimal  fraction  of  the  observation 
noise  power  ( cfo)  and  I  is  the  identity  matrix  i.e.,  ones  only  on  the  principal 
diagonal.  Unfortunately,  one  does  not  know  how  much  noise  to  remove.  Tranter 
[9]  has  studied  this  technique  for  the  cases  of  one  or  two  real-valued  noisy 
sinusoids.  The  Pisarenko  Harmonic  Decomposition  (PHD)  technique  is  a  special 
case  which  will  be  examined  subsequently  in  Section  IV.  It  does  not  assume 
that  the  sinusoids  are  harmonically  related  and  is  a  special  ARMA  process 
which  has  equal  AR  and  moving-average  parameters. 

I.  Computing  R~l  from  LPF  Coefficients 

Although  the  PSD  of  the  AR(M)  process  can  be  determined  directly  from 
the  LPF  coefficient  set  (aM(k)},  i.e.,  without  ever  generating  the  actual 
correlation  matrix  R;  there  are  practical  reasons  for  needing  R-  inverse 
(R~l).  Burg  addressed  this  point  in  his  1968  paper  [12]  and  Andersen  [7]  gave 
a  general  solution  in  terms  of  the  PEF  sets  and  prediction  errors  (Pm)  for 
orders  m  -  0  M.  Burg  showed  in  his  dissertation  [21]  that  the  inverse 
correlation  matrix  could  be  expressed  as 


R-l 


t 

=  ADA 


(24) 
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Direct  solution  of  Eqn.  (20)  by  matrix  inversion  has  been  documented  by 
Barrodaie  and  Erickson  [14]  and  will  be  referred  to  as  the  B|rrodale-Erickson 
Algorithm  (BEA).  Unfortunately,  such  inversion  requires  0(m  )  operations,  and 
Marple  [6]  has  proposed  an2alternative  algorithm  which  exploits  the  structure 
[rm(i,j)]  and  requires  0(m  )  operations.  This  algorithm  will  hereafter  be 
referred  to  as  Marple's  Least-Squares  Algorithm  (MLSA).  The  BEA  and  MLSA 
techniques  eliminate  SLS  and  significantly  reduce  the  bias  and  variance  of  the 
resulting  spectral  estimates.  Additional  observations  regarding  these  two 
least-squares  algorithms  are  included  in  Section  IV. 


F.  AR  Parameter  Estimation  for  Long  Data  Records 

The  Maximum  Likelihood  Estimator  (MLE)  is  essentially  equivalent  for 
N  »  m  to  solving  the  YWE  using  Eqn.  (13)  to  find  estimates  of  R(k). 
Unfortunately,  finding  the  exact  MLE  solution  of  the  AR(m)  parameter  set  is 
difficult.  Consequently,  three  sequential  estimation  schemes  for  updating 
(am(k) }  of  a  slowly  varying  AR(m)  process  are  documented  by  Kay  &  Marple  [2, 
pp.  1393-5).  They  include  a  recursive  least  squares  method  which  resembles 
Kalman  filtering,  an  adaptive  LPF  which  uses  a  gradient  technique  with 
corresponding  slow  covergence,  and  a  lattice  filter  which  provides  updates 
only  for  the  reflection  coefficients  coupled  with  the  LDE  to  update  the 
remaining  m-2  coefficients.  The  applications  contemplated  in  this  study  are 
restricted  to  short  records. 

G.  AR  Model  Order  Selection 

The  technique  employed  in  either  LSSE  or  MESE  has  been  to  increase 
the  LPF  order  (m)  until  the  resulting  error  energy  (E^,)  reaches  some  accep¬ 
table  level.  Unfortunately,  it  follows  from  Eqn.  (11)  that  E^j  <  E„,  and  it 
has  been  shown  for  an  AR(M)  process  that  the  estimated  spectrum  is  highly 
smoothed  when  m  «  M  and  spurious  detail  results  when  m  »  M.  Akaike  has  pro¬ 
posed  two  prediction  estimators  which  have  received  considerable  atten¬ 
tion,  the  Final  Prediction  Error  (FPE)  and  the  Akaike  Information  Criterion 
(AIC)  which  for  a  mc^  -  order  LPF  are  given  by 

and  (22) 

AlCn,  -  ln(Em)  +  2(m+l)/N  . 

A  third  estimator  due  to  Parzen  [15]  called  the  Criterion  Autoregressive 
Transfer  (CAT)  is  given  by 


errors  (presumed  to  be  known).  The  reflection  coefficients  so  derived  are 
given  by 


f 


Ktn+1 


N-l  * 

E  bm  (i-1)  fm(i) 

_2  i*m+l _ 

N-l 

E  {  |  bm(i-l)  |  +  |fm(i)  |2  } 

i=m+l 


(17) 


where  Andersen  [7]  provided  a  flowchart  for  the  Burg  algorithm,  and  in  a  later 
paper  [13]  provided  a  recursive  relation  for  the  denominator  (D^j)  of  Eqn. 
(17), 

Dm+1  -  Dm.(l-  K  |2)  -  |bn(N-m-l)  |2  -  |  fm(m+l)  |2  (18) 

The  computer  program  based  on  the  original  Andersen  flowchart  [7]  was  used  by 
Tranter  [8]  and  will  hereafter  be  referred  to  as  Andersen's  Burg  Algorithm 
(ABA).  The  ABA  technique  has  several  problems  including  spectral  line 
splitting  (SLS)  and  significant  bias  and  variance  in  the  PSD  estimate. 


E.  Least-Squares  Spectral  Estimation  (LSSE) 


The  LSSE  technique  removes  the  LDE  constraint  and  partially  differen¬ 
tiates  Eqn.  (16)  with  respect  to  all  the  (atn(k)},  setting  the  results  equal  to 
zero,  which  yields 

m 

2  E  ajj)  rm(i,j)  „  0 
aa^Ci)  j-0 

where  (19) 

N-m-1 

rm(i,1)  -  E  x(k+m-j)  x*  (k+m-i) 
k-0 

+  x(k+i)  x*  (k+j ) . 

In  matrix  form  Eqn.  (19)  becomes 

t  rm^  )  J(m+1  ,m+l)  l  ®m^^  (m+1 ,1)  (20) 

”  ^Em  ®  (1  ,m+l) 

where  the  minimum  prediction  error  energy  is 


m 

Em  -  T  VJ)  rm(0  ,  J  )  ♦ 


(21) 


D.  Forward-Backward  Linear  Prediction  (FBLP) 

The  FBLP  technique  was  proposed  by  Burg  [12]  with  forward  and  back¬ 
ward  PEF  errors  given  by 

m 

fm(n)  =>  x(n)  +  E  a(i)x(n-i) 

i-1  (14) 

m 

bm(n)  =  x(n-m)  +  E  a*(i)x(n-m+i)  , 

i-1 

where  fm(n)  =  x(n)  -  x(n) ,  x(n)  given  by  Eqn.  (2)  with  all  (b(i)}  *  0. 

Likewise  bm(n) ,  no  relation  to  b(i),  represents  a  backward  prediction  error 
estimate  of  x(n-m)  based  on  "future"  samples.  Note  that  the  summation  limits 
in  Eqn.  (14)  are  restricted  such  that  no  samples  of  x(n)  outside  the  range  n»0 
to  N  are  involved,  i.e.  ,  no  assumptions  are  made  regarding  the  time  series 
outside  the  N-sample  interval;  in  contrast  to  Eqn.  (13)  where  R(k)  is  esti¬ 
mated  based  on  the  tacit  assumption  that  x(n)  =  0  for  n  <  0  or  n  >  N.  The 
forward-backward  errors  can  be  expressed  by  the  recursion  relations 


fm<n>  *  fm-l<n)  +  W-l*®”1) 

* 

bn/n>  *  bm-l(n-D  +  V®-l(*0 


(15) 


subject  to  the  initial  conditions  fQ(n)  *  b0(n)  »  x(n).  Implemented  in  a 
flowgraph  form  with  error  outputs  from  the  previous  section  as  inputs  forms  a 
so-called  lattice  PEF,  two  sections  of  which  are  shown  in  Figure  1  where  (•) 
implies  a  summing  node.  To  obtain  estimates  for  {aj^(k) >  Burg  minimized  the 
total  error  energy 


Z  l  K<«>  |2  +  K<“>  |21 


subject  to  the  LDE  constraint  for  a(i),  thus  insuring  a  stable  AR  filter. 
Setting  the  derivative  of  Em,  m  =  1,2,  . ..,  M,  with  respect  to  K^j  equal  to 
zero  permits  solving  for  {Km}  as  functions  of  the  prior-order  prediction 


Figure  1.  Lattice  PEF. 


Af 


1 


cos(2<t>  +  BN)  sin  (ON) 


(43) 


2  ttNT 

It  is  apparent  that  the  error  is  a  decreasing  function  of  the  number  of 
samples  used;  however,  when  normalized  with  respect  to  the  nominal  DFT  resolu¬ 
tion  celi  width  (1/NT),  the  maximum  error  is  roughly  16%  of  this  cell  width 
and  occurs  for  four  values  of  6.  Swingler  [46]  showed  that  the  frequency  error 
could  be  reduced  by  inserting  into  the  reflection  coefficient  expression,  Eqn. 
(17),  a  Hamming  weight  term  in  both  numerator  and  denominator  summations. 
Moreover,  independently  of  both  Nuttall  [47]  and  Barrodale  &  Erickson  [14]  he 
also  proposed  the  least  squares  approach  which  is  implemented  by  Marple's 
algorithm.  A  good  review  of  noisy  sinusoid  spectral  estimation  sensitivity  of 
the  ABA  to  such  factors  as  sample  size,  order  number,  phase,  and  SNR  is 
detailed  by  Chen  &  Stegen  [48].  It  has  been  shown  by  Marple  [6]  that  the  MLSA 
virtually  eliminates  both  frequency  bias  as  a  function  of  <t>  and  SLS  associated 
with  large  SNR.  However,  little  additional  has  been  published  regarding 
sample  size  effects.  The  selection  of  order  is  handled  routinely  in  MSLA  by 
specifying  MM,  after  which  the  Marple  algorithm  increases  M  until  either  one 
of  the  tolerance  tests  on  error  energy  ratios  is  met  or  M=MM.  The  Akaike  FPE 
first  minimum  test  is  incorporated  into  the  BEA  [14]  and  Tranter  [8]  evidently 
used  it  in  deciding  the  order  for  ABA  too. 

Landers  &  Lacoss  [17]  noted  that  the  various  prediction  error 
algorithms  were  somewhat  sensitive  to  real  vs.  complex  data  in  the  presence  of 
uniform  white  noise.  The  results  for  a  noisy  sinusoid  using  the  ABA  and  MLSA 
are  summarized  for  one  particular  test  in  Table  1.  The  value  of  M  selected  by 
the  MLSA  and  corresponding  frequency  estimate  (F  =  fcT  =  0.05)  for  32  samples 
of  a  noisy  sinusoid  with  MM  =  8  are  tabulated  for  both  real  and  complex  data. 
The  test  was  repeated  for  the  ABA  using  only  complex  data  (since  Af  would  be  a 
maximum  at  <b  =  50°),  initially  with  M=8  and  then  M=3  which  was  the  order  pre¬ 
dicted  by  the  first-minimum  of  both  the  FPE  and  CAT  criteria.  (The  AIC  had 
not  reached  a  first  minimum  by  M=8.)  To  illustrate  the  effect  of  changing  N, 
the  10  dB  SNR  case  was  run  for  complex  data  with  N  reduced  from  32  to  16.  The 
spectral  estimates  (M=8)  for  the  ABA  and  MLSA  were  0.0500  and  0.0506,  respec¬ 
tively.  Doubling  N  from  32  to  64  yielded  0.0488  and  0.0497,  respectively. 
These  results  would  tend  to  show  that  a  noisy  sinusoid  can  be  reasonably  esti¬ 
mated  for  SNR  >  0  dB  even  with  a  small  number  of  samples.  Comparative  results 
from  ABA  and  MLSA  for  real-valued  noisy  sinusoid  (fcT  =  0.0725,  a  =  50,  d>  = 
50°)  are  shown  in  Figures  A-3  and  A-4.  The  ABA  frequency  bias  for  this  par¬ 
ticular  <t>  is  NTAf  =  0.125  (12.5%).  It  is  evident  from  Table  1  that  the  ini¬ 
tial  phase  is  not  a  problem  with  complex-valued  samples. 
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Table  1.  MLSA  and  ABA  Performance  Comparison 
Noisy  Sinusoid  (fcT  *  0.05,  A  -  50°,  N  *  32) 


SNR 

MLSA  (MM  = 

Real 

M  F 

8) 

Complex 

M  F 

ABA 

Complex 

F(M=8)  F(M=3) 

0. 

8 

0.0494 

8 

0.0455 

0.0476 

0.0440 

10. 

8 

0.0494 

8 

0.0485 

0.0497 

0.0494 

20. 

8 

0.4097 

8 

0.0497 

0.0500 

0.0500 

30. 

6 

0.0500 

4 

0.0500 

0.0500 

0.0500 

40. 

2 

0.0500 

2 

0.0500 

0.0500 

0.0500 

C.  Closely-Spaced-Sinusoids  Estimation 

The  ability  to  resolve  two  closely-spaced  noisy  sinusoids  has 
received  considerable  attention  in  the  literature.  Marple  [49]  empirically 
determined  that  the  ability  to  resolve  two  sinusoids  at  frequencies  fj  and 
f2  depended  on  the  PEF  order  (M)  and  SNR  (a,  not  in  dB) ,  i.e., 

i  i  1.023 

^  -  I  f2  “f 1  I  T  *  M[ a(M+l) jO-31  '  (44) 

This  formula  was  obtained  assuming  M  correlation  lags  and  is  reproduced  in 
Table  2  for  M  ranging  from  1  to  30  and  SNR  ranging  from  0  to  50  dB.  If  the 
number  of  samples  N  »  M,  Eqn.  (44)  serves  as  an  upper  bound  on  resolution, 
i.e,  increasing  M  by  two  or  three  should  yield  equivalen  resolution  for  a 
given  SNR.  Marple  also  observed  that  a  signal  of  the  form 

x(k)  =  Aj  cos(0]k  +  <p j )  +  A2  cos(02k  +  *2)  (45) 


was  most  difficult  to  resolve  when  <p ^  =»  d>2.  The  frequency  resolution  when 
compared  to  the  conventional  DFT  periodogram  was  roughly  four  times  as  good  as 
20  dB,  twice  at  0  dB ,  and  equivalent  at  -10  dB. 

Barrodale  &  Erickson  [14]  compared  their  BEA  with  the  conventional 
ABA  and  showed  that  it  outperformed  the  ABA  for  a  variety  of  M  choices  and  a 
small  number  of  samples.  Specifically,  for  N=75  and  fj  =  0.003,  f2  *  0.02 
(which  corresponds  to  0.23  and  1.5  cycles,  respectively)  the  BEA  program  was 
able  to  identify  both  f[  and  f2  correctly  for  M=16  and  f2  for  M  5,  whereas 
the  ABA  generated  false  spectral  estimates  for  all  examples.  However,  they 
observed  comparable  performance  if  at  least  15  cycles  of  sinusoid  samples  were 
available.  Ulrych  &  Clayton  [50]  showed  that  their  least-squares  algorithm, 
the  basis  for  both  BEA  and  MLSA,  had  considerably  less  frequency  bias  than  the 
ABA.  They  also  recommended  that  M  be  selected  in  the  range  N/3  <  M  <  N/2 
rather  than  relying  on  the  FPE  or  AIC  first-minimum  concept.  This  recommen- 


Table  2.  Frequency  Resolution  Estimates  Using  Equation  44. 


s m  iD»> 


N 

t. 

14. 

24. 

34. 

44. 

54. 

1 

4.93*75 

4.44*99 

4.19924 

9.497*1 

4.44744 

4.42341 

t 

t. 3*431 

4.17941 

4.44747 

4.94344 

4.42144 

4.91932 

3 

t. 21337 

4.14944 

9.4S3S4 

4.92*24 

9.41245 

9.44*34 

4 

t. 19433 

4.47*97 

9.43794 

4.41437 

4.44944 

9.4*441 

9 

t.lltlt 

4.4S749 

9.42439 

4.41349 

4.44*49 

4.44333 

• 

t .49399 

4.44999 

4.42252 

4.41143 

4.44544 

(.442*5 

7 

4.97722 

4.43792 

4.41452 

4.44947 

9.44444 

4.44214 

• 

t.ttsis 

4.43191 

4.415*3 

4.447*5 

4.44375 

4.44144 

9 

t.tcetc 

4.42745 

4.41344 

4. 44*91 

4.44323 

9.44154 

19 

4.44997 

4.42399 

4.41179 

4.4*575 

4.44242 

4.4*134 

11 

9.94334 

4.42123 

4.41444 

4.495*9 

4.44249 

9.44122 

It 

4.43t79 

4.4199t 

4.49934 

4.44495 

4.4*223 

4.44199 

13 

4.4349* 

4.91712 

4.49439 

4.94411 

4.4*241 

9.44499 

14 

9.93177 

4.41SS4 

4.447*2 

4.44373 

4.44193 

4.44494 

1* 

4.92997 

4.41424 

4.44*97 

4.44342 

4.441(7 

4.44442 

It 

4.92*74 

4.41314 

4.99*42 

4.44314 

4.4*154 

9.44*75 

17 

9.92473 

4.41211 

4.49993 

4.44291 

4.49142 

4.44*74 

It 

4.92297 

4.41125 

4.49951 

4.4*274 

9.44132 

4.444*5 

19 

4.92142 

4.41449 

4.4*914 

9.4*252 

4.44123 

(.4*4*4 

at 

4.92*94 

4.44941 

4.4*441 

4.44239 

4.44115 

(.4*456 

ti 

4.41991 

4.44921 

4.4*451 

9.4*221 

9.94144 

(.44453 

at 

4.41771 

9.449*7 

4.49429 

4.4*244 

4.44142 

(.44*54 

13 

4.91*72 

9.44919 

4.4*441 

4.4*199 

9.4449* 

*.***•7 

14 

4.41991 

4.44779 

9.44394 

4.9914* 

9.4*491 

9.4*445 

at 

4.41S44 

4.94739 

4.443*4 

4.9*179 

4. 44*4* 

4.44*42 

at 

4.91424 

4.44*94 

9.44342 

9.4*164 

9.4*442 

*.***44 

t? 

4.91391 

4.94M5 

4.4432* 

4.441*4 

9.44474 

4.44*34 

aa 

4.41199 

9.4404 

4.44311 

9.44152 

4.4*475 

4.4*43* 

aa 

4.4107 

9.44*94 

4.44297 

9.44145 

4.44471 

4.44435 

at 

9.91194 

4.44944 

4.9429* 

9.94139 

9.944*4 

4.44*33 

21 


dation  is  also  supported  by  Kay  &  MarpLe  [2,  p.  1397]  who  illustrated  that 
increasing  M  from  four  to  32  resolved  two  sinusoids  at  0.143  and  0.200  at  a 
SNR  3  0  dB.  An  example  [16]  of  spectral  estimation  for  two  equal  amplitude, 
complex-valued  sinusoids  (fj  3  0.11,  f2  =  0.14,  N=40,  a=50)  using  the  ABA 
program  with  M=15  is  given  in  Figure  5.  Spectral  estimates  for  the  two  peaks 
differ  slightly,  e.g.  ,  ABA:  (0.1120,  0.1388)  and  MLSA:  (0.1105,  0.1395). 

The  PEF-order  estimators  for  the  ABA  program  are  shown  in  Table  3.  Local 
minima  for  the  FPE  and  CAT  estimators  track  each  other  rather  consistently, 
whereas  the  AIC  estimate  does  not  provide  a  local  minimum  through  M=*15,  which 
is  consistent  with  the  MLSA  setting  M  =  MM  =  15. 

Table  3.  PEF  Order  Estimators  Using  ABA  Program  (M=15). 

Two  Noisy  Sinusoids  (fjT  3  0.11,  f2T  3  0.14,  SNR  =  17  dB) 


PCUER  -  PU>  A<A1KE  -  fPE(J)  AroUE  -  Pi:..)  e~;;'LN  -  ■ 

0.96O4051E-01  0.10s:  sme  •**  -9  6  37:9«ie*pc 
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0.3I02S25E-01  0.7239226E-01  -0  1085U2E  »93  1433034>E-«>; 


Tranter  [8]  developed  a  real-valued  (N=100)  test  case  for  F.qn.  (45) 
in  which  the  chosen  frequencies  (fj  *  0.098,  f2  •  0.102)  were  spaced  less  than 
one-half  a  DFT  cell  apart.  Both  he  and  Shumway  [10]  used  the  BEA  program  in 
conjunction  with  the  FPE  first-minimum  criterion  to  study  the  effects  of  SNR 
on  the  spectral  estimation  and  concluded  that  proper  resolution  required  at 
least  40  dB  SNR.  This  value  turns  out  to  be  too  large  because  of  their 
reliance  on  the  FPE  first-minimum  criterion  to  select  M. 

The  results  obtained  by  repeating  this  test  case  using  MLSA  on  real 
data  for  a  variety  of  SNR-M  combinations  are  summarized  in  Table  4  and  some 
responses  are  shown  in  Figure  6.  It  is  apparent  that  the  test  frequencies  can 
be  resolved  with  less  than  1%  error  for  a  SNR  of  20  dB  provided  M  is  large 
enough.  In  contrast,  an  M  of  7  can  be  used  when  the  SNR  3  50  dB.  The  third 
column  of  Table  4  labeled  S(0.1)  indicates  the  response  in  dB  at  the  mid-point 
(F=»0.1)  relative  to  the  smaller  sinusoid,  with  0  dB  indicating  failure  to 
resolve  the  two  components.  The  final  column  (Me)  indicates  the  value  of  M 
predicted  from  Table  2  which  should  yield  a  response  similar  to  the  33/20  dB 
curve  in  Figure  6  for  that  particular  SNR.  As  stated  earlier  Me  is  based  upon 
autocorrelation  estimates  rather  than  data  samples;  however,  the  difference 
between  it  and  the  M  which  yields  S(0.1)  3  3  dB  is  evidently  a  function  of 
SNR,  ranging  from  approximately  nine  at  20  dB  to  three  at  50  dB. 
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Figure  5.  Double  sinusoid  spectral  estimate  using  ABA 
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Table  4.  Summary  of  Tranter's  Test  Case  Using  MLSA 


M 

SNR  (dB) 

S(0.1)  dB 

Me 

33 

20 

2.5 

24 

33 

30 

24.0 

14 

16 

30 

0. 

14 

16 

40 

10.5 

8 

10 

50 

18.5 

5 

7 

50 

2.0 

5 

5 

50 

0. 

5 

D.  Target  Identification 

Target  Identification  is  a  potential  radar  application  for  the  abi¬ 
lity  of  MSA  techniques  to  resolve  closely-spaced  sinusoids.  Unfortunately, 
little  has  been  found  In  the  open  literature  regarding  this  feature  and  time 
did  not  permit  a  search  of  the  classified  literature.  Taylor  £t.  al.  [51] 
described  one  situation  in  which  the  DFT  and  MEM  techniques  were  both  applied 
to  two  sweep-scan  radar-return  data  sets  (N-64)  from  a  single-engine  jet 
aircraft.  Basically,  the  MEM  spectrum  was  able  to  identify  spectral  lines  in 
both  data  sets  for  the  target  doppler  and  upper  modulation  sideband  (caused  by 
engine  turbine  rotation),  whereas  the  DFT  spectrum  was  poorer  and  failed  to 
identify  the  modulation  sideband  for  one  data  set.  A  classified  report  by 
Gardner  [52]  indicates  that  more  complicated  spectral  returns  would  be  needed 
to  distinquish  various  types  of  aircraft.  Only  Swingler  [45]  seems  to  have 
addressed  the  resolution  of  more  than  two  sinusoids,  providing  some  infor¬ 
mation  regarding  the  resolvability  of  ten  harmonically  related  sinusoids  with 
relative  amplitudes  in  the  range  0.1  to  1.0  with  5%  added  white  noise.  Time 
did  not  permit  simulation  of  this  situation  nor  the  review  of  some  specialized 
techniques  to  resolve  sinusoids  in  the  presence  of  noise.  These  methods, 
reviewed  by  Kay  &  Marple  [2,  pp.  1402-7]  include  Pisarenko  Harmonic 
Decomposition,  Prony  Spectral  Line  Decomposition  and  Prony’s  Method  (Extended) 
with  appropriate  references  cited.  In  each  case  the  number  of  sinusoids  must 
be  estimated  to  establish  a  model  order.  Other  recent  papers  on  this  topic 
have  been  written  by  Tufts  &  Kumaresan  [53,  54]  Kay  [55],  Marple  [56,  57],  and 
Shumway  [ 10] . 


V.  RESULTS  AND  RECOMMENDATIONS 


n 


The  basic  theory  of  MSA  was  traced  from  the  Yule-Walker  equations,  which 
assumed  the  autocorrelation  function  was  either  known  or  could  be  approxi¬ 
mated,  thru  the  forward-backward ' energy  minimization  of  Burg,  to  the 
algorithms  of  Andersen,  Barrodale  &  Ericksen,  and  Marple  which  are  used  in 
MSA.  The  various  model-order  selection  schemes  were  used  with  both  AR  and 
noisy-sinusoid  signals.  Sawyer's  method  was  recorded  for  determining  the 
inverse  of  the  correlation  matrix  using  iterative  equations  which  are  func¬ 
tions  of  the  Mc^  order  LPF  coefficients. 

Clutter  modeling  and  subsequent  suppression  by  inverse  filtering  seems  to 
be  a  definite  application  of  MSA.  It  was  shown  that  clutter  is  basically 
modeled  spectrally  as  an  AR  process  and  the  McMaster  University  team  (Haykin, 
et.  al. )  demonstrated  that  both  simulated  and  real  clutter  PSD  could  be 
modeled  adequately  by  MSA,  particularly  when  the  samples  were  complex  valued. 
Using  theory  developed  in  the  late  1960s,  Sawyers  demonstrated  the  feasibility 
of  designing  a  DFB,  which  would  optimally  improve  the  SIR,  by  multiplying  the 
inverse-R  matrix  by  a  steering  vector.  This  matrix  could  be  obtained  by  tra- 
diitonal  means,  if  the  clutter  correlation  function  was  specified,  or  by  the 
use  of  MSA  as  discussed  earlier  when  only  clutter-dominated  sample  returns 
were  available.  Adaptive  filtering  based  on  this  technique  has  been  studied 
also  by  Nitzberg  in  support  of  a  U.  S.  Air  Force  RADC  project.  Gibson  & 

Haykin  have  studied  properties  of  adaptive  lattice  PEF  required  to  suppress 
clutter.  Such  filters  require  computation  only  of  reflection  coefficients 
obtained  from  Burg's  forward  and  backward  errors.  Decoupling  appears  to  be  a 
problem  since  all  coefficients  initially  adapt  to  the  strongest  spectral  com¬ 
ponent  and  once  Kj  is  fixed,  the  remaining  values  readjust  to  the  next 
strongest  component,  etc.  Their  research  indicated  roughly  three  stages  for 
every  broad  peak  in  the  PSD  and  one  for  each  sharp  peak. 

Resolution  of  two  sinusoids  was  studied  for  a  variety  of  SNR  values.  It 
was  shown  that  the  conventional  model-order  predictors  of  Akaike  and  Parzen 
tended  to  produce  a  first-minimum  at  too  low  an  order.  Marple' s  technique 
seems  to  be  more  appropriate;  however,  the  tolerance  choices  must  be  selected 
with  some  care.  Marple' s  empirical  formula  for  frequency  resolution  capa¬ 
bility  as  a  function  of  model  order  and  SNR,  as  reproduced  in  Table  2,  seems 
to  provide  a  reasonable  starting  point  for  selecting  M  for  a  particular  SNR. 
Studies  performed  in  Section  IV  tend  to  indicate  that  true  AR(M)  processes  can 
be  spectrally  estimated  quite  well  with  m  =  M  and  noisy  sinusoid  pairs  can  be 
resolved  with  large  m,  typically  approaching  N/3  for  low  SNR. 

It  is  recommended  that  clutter  suppression  algorithm  be  designed  with 
Ouiet  Radar  parameter  sets  and  tested  with  complex  data  obtained  from  the 
radar.  It  is  also  suggested  that  the  MLSA  be  applied  to  the  problem  of 
multiple  sinusoid  estimation  where  the  amplitude-frequency  assignment  is  cho¬ 
sen  to  emulate  the  spectrum  associated  with  a  particular  type  or  class  of 
targets.  Finally,  contacts  should  be  Initiated  with  researchers  at  McMaster 
University,  RADC,  and  in  other  organizations  where  work  is  underway  on  the 
clutter  suppression  problem. 
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ACRONYM  LIST 


ABA  Andersen's  Burg  Algorithm 

AIC  Akaike  Information  Criterion 

AR  Autoregressive 

AR(m)  AR  process  of  mtll-order 

ARMA  Autoregressive  Moving  Average 

BEA  Barrodale-Erickson  Algorithm 

CNR  Ciutter-to-Noise  Ratio 

DFB  Doppler  Filter  Bank 

OFT  Discrete  Fourier  Transform 

FBLP  Forward-Backward  Linear  Prediction 

FPE  Final  Prediction  Error 

LDE  Levinson-Durbin  Equation 

LPF  Linear  Prediction  Filter 

L.SSE  Least-Squares  Spectral  Estimation 

MEM  Maximum  Entropy  Method 

MESE  Maximum  Entropy  Spectral  Estimation 

MLE  Maxirauim  Likelihood  Estimator 

MLSA  Marple's  Least-Squares  Algorithm 

MSA  Modern  Spectral  Analysis 

MTI  Moving  Target  Indicator 

PDF  Probability  Density  Function 

PEF  Prediction  Error  Filter 

PHD  Pisarenko  Harmonic  Decomposition 

PSD  Power  Spectral  Density 

SCR  Signal-to-Clutter  Ratio 
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ACRONYM  LIST  (Cont’d) 


Signal-to-Int4rf erence  Ratio 
Spectral  Line  Splitting' 
Signal-to-Noise  Ratio 
Yule-Walker  Equations 


33(34  Blank) 


APPENDIX 


MSA  PROGRAM 


The  computer  progam  listed  in  Figure  A-l  of  this  appendix  is  basically  an 
extension  to  complex  data  of  Tranter's  progam  [8].  The  complex-data  extension 
is  based  on  the  algorithm  of  Haykin  &  Kesler  [58]  ,  Datta  [59]  and  Bowling 
[60].  The  BEA  [14]  used  by  Tranter  for  LSSE  was  replaced  by  Marple's 
algorithm  [6]  because  it  is  computationally  more  efficient  and  has  its  own 
prediction-order  estimator,  which  was  shown  in  Section  IV  to  be  better  than 
the  conventional  estimators  of  Akaike  and  Parzen  used  herein  in  conjunction 
with  the  ABA  option.  Marple  estimates  that  his  LSSE  algorithm  requires  only 
20%  more  computation  time  than  the  ABA  program,  whereas  Tranter  estimated  that 
the  BEA  program  took  three  times  as  long.  An  X-Y  plot  routine  written  by  Tom 
Cash  for  the  PDP  11/34  computer  was  also  included  in  the  MSA  program. 

The  main  program  is  designed  to  read  user  control  parameters  (see  Table 
A-l),  generate  the  real  (imaginary)  noise  power  a  consistent  with  the  user 
supplied  SNR  and  Eqn.  (42),  call  the  appropriate  subroutines  (S.R.),  remove 
any  bias  (mean  value)  from  the  generated  samples,  write  appropriate  output 
information  information,  and  generate  PSD  (dB)  vs.  frequency  (Hz)  data  points 
for  S.R.  PLT.  The  various  subroutines  and  their  purposes  are  itemized  in 
Table  A-l  in  the  order  in  which  they  are  listed  in  Figure  A-l. 

Some  examples  of  the  printouts  and  PSD  plots  achieved  with  the  ABA  and 
MLSA  follow.  The  common  header  describing  the  test  case  data  is  deleted  from 
the  b-part  of  Figures  A-2  and  A-3  to  insure  a  one-page  fit. 


Table 

A-l. 

Input  Parameters  for  Program  MSA. 

Card 

Symbols 

Function 

1 

ISIG 

Test  flag  indicating  presence  (»1)  or  absence 
(-0)  of  sinusoids,  Eqn.  (45). 

ICLT 

Ibid,  for  AR(MC),  Eqn.  (2). 

ICPX 

Data  samples  are  real  (»0)  or  complex  (*1). 

IALG 

Estimation  technique  is  ABA  (=0)  or  MSLA  (=*1). 

2 

SNRDB 

SNR  (dB)  ( ISIG=1 ) . 

CNRDB 

CNR  (dB)  (iCLT=l) . 

3 

Al,  FI, 

PHI1 

Amplitude,  frequency,  and  initial  phase  of 
first  sinusoid  in  Eqn.  (45).  (ISIG=1). 

4 

A2 ,  F2 , 

PHI2 

Ibid,  for  second  sinusoid. 

5 

N 

Number  of  samples  in  test  waveform. 

A-3 


Figure  A-l .  Computer  program  for  MSA.  (Cont'd) 


Figure  A-l .  Computer  program  for  MSA.  (Cont'd) 
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B)  ABA 


Figure  A-3 . 


Example  printouts  for  noisy  sinusoid. 


par Zen  -  cot  : > 


-*.*3349180*1 
-9.2543638092 
*9. 4353422(482 
-9.541839X482 
-9.643885X.92 
-9.6766258E.92 
-9. 68 13776092 
-9. 6672352E»82 
-9.6541871E482 
-9.6589137(492 
-9. 643698%  *92 
-8.6322446E.82 
-9.6441886C.92 
-9. 629889 4E*82 
-8.6157121E.82 
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Table  A-l.  Input  Parameters  for  Program  MSA.  (Cont'd) 


Card 


Symbols 


Function 


6 


7 


M 

MC 


Maximum  order  of  AR  Prediction  used  to  ter¬ 
minate  ABA  and  possibly  MLSA. 

Order  of  AR  test  waveform. 


K 

FN 

FI,  FF 

TOL  1 


TOL  2 


AR 


Number  of  points  in  PSD  plot. 

Sampling  frequency  (Hz) . 

Initial  and  final  frequencies  in  PSD 
plot  (Hz). 

Ratio  of  final  to  initial  estimation  error 
in  MLSA  which  if  reached  will  terminate 
search.  (Default  -  0.001). 

Ratio  of  present  estimation  error  to  previous 
order  error  in  MLSA  which  if  reached  will 
terminate  search.  (Default  =  0.005). 

Array  of  MC  coefficients  for  AR(MC)  test 
signal  (ICLT  -  1). _ 


Table  A-2.  Subroutines  for  MSA  Program. 


S.  R. 

SIG1 

GASS 

NPWR 

FPE1 

MFREQR 

PLT 

MEMC 

ARCLTR 

MARPLE 


Purpose 

Generate  sinusoidal  samples,  Eqn.  (45). 

Generate  Gaussian  noise  samples  with  zero  mean 
and  standard  deviation  o. 

Calculate  sample  mean-square-value  for  signal, 
noise,  and  signal  plus  noise. 

Calculate  three  prediction-error  estimates 
( IALG  *  0),  Eqns .  (22)  and  (23). 

Compute  PSD,  Eqn.  (5). 

Generate  plot  of  PSD  (dB)  vs.  Frequency  (Hz) 
and  identify  frequency  location  of  spectal  peak. 

Complex-data  version  of  ABA.  (IALG  ■  0) . 

Generate  AR(MC)  samples,  Eqn.  (2). 

Marple's  least-squares  algorithm.  (IALG  *  1). 
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Example  1  =  AR(4)  Model  (CNR  =  50dB) 

This  example  estimates  the  PSD  for  a  given  AR(4)  process  with  coefficients 
shown  in  “AR  CLUTTER  COEFF.  ="  line  in  Figure  A-2(a).  The  Marple  algorithm 
(MLSA)  terminated  with  M=4  (the  proper  choice)  because  E(4)/E(  <+>)  <  0.00]  as 
shown.  The  AR  coefficient  closely  approximate  the  negatives  of  the  defined 
values  because  the  minus  sign  of  Eqn.  (2)  was  inadvertently  omitted  from  S.R. 
ARCLTR.  This  could  readily  be  corrected  by  the  change 

XC(J)  =  CMP  LX  (-X(K),  $.<]>). 

The  ABA  was  permitted  to  run  to  M=8  to  show  the  three  predictor-order 
results.  Using  first  minimum  choices,  the  proper  order  would  be  FPE=5,  AIC 
8  (with  a  rather  flat  value  for  5  £  M  <_  7),  and  CAT=5.  The  PSD  plot  appears 
in  Figure  A-4. 

Example  2  -  Noisy  Sinusoid  (SNR  =  17  dB) 

The  second  example  is  for  a  noisy  sinusoid  described  by  Eqn.  (41).  The 
user  supplies  the  SNR  (17  dB)  which  is  converted  to  a  number  from  which 
a  (SIGMA)  can  be  determined  in  accord  with  Eqn.  (42).  Note  that  the  second 
signal  amplitude  (A2)  must  be  set  equal  to  zero  as  shown  in  Figure  A-3(a).  In 
this  case  the  prediction  order  reaches  the  preset  maximum  (MMAX)  of  15.  The 
ABA  predictors  shown  in  Figure  A-3(b)  indicate  FPE»7,  AIC  15  (with  flat 
values  in  the  ranges  7  M  9  and  13  ^  M  15),  and  CAT*7.  The  spectral  fre¬ 
quency  for  ABA  is  7.375  Hz  vs.  the  correct  MLSA  estimate  of  7.250  Hz.  The  ABA 
spectral  estimate  shown  in  Figure  A-4  also  has  considerably  more  sidelobe 
energy  over  the  region  dc  to  1/4T. 

Generating  data  samples  for  a  pair  of  sinusoids  simply  requires  supplying 
non-zero  values  for  the  three  parameters  of  Card  4  in  Table  A-l.  The  user- 
supplied  SNR  is  based  on  the  formula 

a  -  (A2  +  A2)/2 a2 
1  2 

and  is  used  for  both  real  and  complex  valued  signals.  Results  of  such  simula¬ 
tions  are  recorded  in  Section  IV.  C. 
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